Report
Gli statistichi?
14 giugno 2018
Indice
2.2. Linear Discriminant Analysis
2.3. Quadratic Discriminant Analysis
1. Introduzione
Lo scopo di questo report è portare avanti un’analisi di classificazione a partire dai dati riguardanti alcune caratteristiche fisico-chimiche dei vini dell’azienda portoghese Vinho Verde. I dati a disposizione riguardando in particolare due tipologie di vino (rosso e bianco) ai quali corrispondono due differenti dataset che verranno uniti in un unico set di dati che comprenderà un totale di 6497 unità statistiche. L’obiettivo dell’analisi è quello di classificare la qualità del vino che nei dati di partenza è rappresentata da una variabile categorica ordinabile con modalità che vanno da 3 a 9.
E’ possibile notare come le modalità 5 e 6 siano le più numerose e che da sole caratterizzino il 76,6% del totale delle unità statistiche, per questo motivo, per ottenere classi bilanciate, viene generata una nuova variabile target a partire dalla variabile Quality definita come segue:
\[ \begin{equation} \textit{binary_quality} = \begin{cases} 1 \textit{se quality} > 5 \\ 0 \textit{altrimenti} \end{cases} \end{equation}\]
Per ogni tipologia di vino sono disponibili le seguenti caratteristiche fisico-chimiche:
- Acidità fissa
- Acidità volatile
- Acido citrico
- Zucchero residuo
- Cloruri
- Anidride solforosa libera
- Anidride solforosa totale
- Densità
- pH
- Solfati
- Alcol
di cui si riportano in forma aggregata le densità:
inserire commenti delle differenze tra rosso e bianco
Sebbene la correlazione non dica nulla sul possibile collegamento causa-effetto delle variabili considerate, può comunque essere utile in ottica esplorativa, indagare la correlazione tra la variabile Quality e i predittori.
Si riporta la matrice di correlazione calcolate considerando l’intero dataset:
La variabile quality (nella sua forma originale) è molto poco correlata con i diversi predittori, solo in 2 casi si ottiene un coefficiente di correlazione lineare maggiore di 0.3 (alcol, density). Alcuni predittori risultano molto correlati tra di loro, in particolare i valori più alti (in valore assoluto) che è possibile osservare nella matrice di correlazione sono:
- Anidride solforosa libera - Anidride solforosa totale (0.72)
- Density-alcol (-0.687)
- density - zuccheri residui (0.553)
Considerando la natura della variabile Quality è opportuno valutare la correlazione tra la variabile Quality e le diverse variabili esplicative considerando il coefficiente di Spearman piuttosto che il coefficiente di Pearson:
Il coefficiente di correlazione di Spearman misura il grado di relazione tra due variabili per le quali non si fa altra ipotesi della misura ordinale, ma possibilmente continua. Diversamente dal coefficiente di correlazione lineare di Pearson, non misura una relazione lineare anche qualora vengano usate misure intervallari.
Per maggiori informazioni sui dati è possibile consultare il sito https://www.kaggle.com/c/uci-wine-quality-dataset
Entrambi gli approcci sono stati condotti sfruttando una divisione dei dati in training (60%) e test (40%) set.
2. Classificazione
In questa sezione sono presentati e messi a confronto vari modelli di classificazione con lo scopo di individuare quello con le migliori performance di previsione fuori campione. Le variabili consideate sono i 12 regressori contenuti nei due file di input più la variabile type derivante dall’unione dei due dataset e che assume valore 1 per i vini bianchi e 0 per quelli rossi. La variabile di output è anch’essa una lavorazione della variabile originaria quality ottenuta assegnando valore 0 ai vini di scarsa qualità (<=5) e 1 altrimenti.
In prima approssimazione abbiamo preso in considerazione la classe dei metodi lineari di classificazione, i.e. modelli che producono una funzione discriminante (o una sua trasformazione monotona) lineare per ciascuna modalità della variabile di output. In alternativa, sono presentati anche Quadratic Discriminant Analysis (QDA), k-Nearest Neighbors (kNN) e Generalized Additive Models (GAM).
2.1. Linear Probabity Model
Il modello di regressione si inserisce nel contesto dei metodi di classificazione che si basano sulla definizione di una funzione descriminante, \(\delta_k(x)\), dalla quale poi viene derivata l’appartenenza di una osservazione ad una particolare classe. In particolare, la soluzione proposta da questa tecnica è quella di scegliere come funzione discriminante \(E[Y|X]\) che, data la natura dicotomica del problema (i.e. una osservazione o appartiene o non appartiene ad una data modalità), equivale alla probabilità di appartenenza ad una generica classe g, \(P(Y=g|X)\).
In generale, immaginando di avere un problema di classificazione multinomiale (G modalità), questo approccio consiste nell’adottare un modello di regressione lineare in cui la variabile risposta Y è una matrice di variabili indicatrici, ciascuna avente valore 1 se l’osservazione appartiene alla relativa modalità della variabile di output e 0 altrimenti. Nel caso specifico di una classificazione binaria, possiamo tuttavia adattare un solo modello al dataset di training ed utilizzarlo per calcolare uno score che servirà poi per classificare ciascuna osservazione in base ad un dato treshold (da determinare). Si noti che lo score calcolato sarà per costruzione un numero reale e, quindi, malgrado lo si utilizzi per modellare la probabilità di appartenenza ad una classe - di qui il nome Linear Probability Model (LPM) -, non è possibile interpretarlo come una probabilità. Tuttavia, è possibile dimostrare che sommando gli score di una singola osservazione ottenuti dai modelli relativi a ciascuna modalità il risultato è 1. Infatti:
\[\hat{Y}_{n x G} = X_{n x (p+1)} * ({X}^{T} X)_{(p+1) x (p+1)}^{-1} * {X}_{(p+1) x n}^{T} * Y_{n x G}\]
Post-moltiplicando poi per un vettore colonna G-dimensionale formato da soli 1, possiamo sommare gli elementi di \(\hat{Y}\):
\[\hat{Y}_{n x G} * 1_{G x 1} = X_{n x (p+1)} * ({X}^{T} X)_{(p+1) x (p+1)}^{-1} * {X}_{(p+1) x n}^{T} * Y_{n x G} * 1_{G x 1}\]
A questo punto, pre-moltiplicando a destra per la matrice identità di ordine \(n\) ottenuta come \((X {X}^{T})^{-1}(X {X}^{T})\), è facile verificare che l’equazione precedente si riduce semplicemente a:
\[(X {X}^{T})^{-1}(X {X}^{T}) * X_{n x (p+1)} * ({X}^{T} X)_{(p+1) x (p+1)}^{-1} * {X}_{(p+1) x n}^{T} * Y_{n x G} * 1_{G x 1}\]
\[(X {X}^{T})^{-1}X [({X}^{T} * X_{n x (p+1)}) * ({X}^{T} X)_{(p+1) x (p+1)}^{-1}] * {X}_{(p+1) x n}^{T} * Y_{n x G} * 1_{G x 1}\]
\[[(X {X}^{T})^{-1}X * {X}_{(p+1) x n}^{T}] * Y_{n x G} * 1_{G x 1}\]
\[\hat{Y}_{n x G} * 1_{G x 1} = Y_{n x G} * 1_{G x 1}\]
Per cui si ricava che la somma degli score predetti è 1 per ciascuna osservazioni dal momento che le funzioni indicatrici (RHS) vale \(\hat{Y_1}+ ... + \hat{Y_G} = 1\) per costruzione.
Per cui, malgrado gli score non siano probabilità, sono accomunate ad esse da questa proprietà.
Passando all’applicazione, abbiamo fittato il modello utilizzando tutti i regressori disponibili:
codice modello
Come anticipato, lo score predetto dal modello non è in generale compreso tra 0 e 1, come si evince dal seguente grafico.
Per poter classificare le osservazioni in una delle due classi è necessario poi fissare un taglio sul valore predetto dal modello che funga da separatore delle due classi. La scelta di questa soglia non è tuttavia così intuitiva come in altri modelli, proprio perchè l’outcome non è vincolato all’intervallo \([0,1]\). Per questa ragione, la scelta del treshold ottimale è stata fatta sulla base di una ricerca a griglia in modo da minimizzare l’errore di classificazione nel da verificare set:
qui va inserita la parte di ottimizzazione treshold
terminare con risulati di accuratezza e ROC
2.2. Linear Discriminant Analysis
L’analisi discriminante lineare (LDA) consiste in un tentativo alternativo di costruire una funzione discriminante. Infatti, riscrivendo \(P(Y=g|X)\) come \(\frac{f_g(x)\pi_g}{\sum_{l=1}^{G}f_l(x)\pi_l}\) si deduce che per stimare la probabilità di appartenenza ad una classe basta conoscere la funzione di densità della classe stessa (e la sua probabilità a priori). Partendo da questa formulazione del problema, l’analisi discriminante lineare consiste nell’assumere che le \(f_g(x)\) siano gaussiane \(p+1\)-dimensionali con matrice di covarianza comune per tutte le classi.
Una volta stimate le \(f_g(x)\) dai dati campionari, è possibile dimostrare che la funzione discriminante associata ad ogni classe risulta essere: \[\delta_g(x)=x^{T}\Sigma^{-1}\mu_g -\frac{1}{2}\mu_g^T\Sigma^{-1}\mu_g+\log{\pi_g}\] Per classificare basta poi scegliere \(g_*=argmax_g\delta_g(x)\)
Risultati e plot distribuzione nella direzione discriminante
Alternativamente, il modello LDA può essere ricavato come combinazione lineare dei regressori che meglio separa le modalità della variabile di output. In altre parole, si cerca \(a^Tx\) tale che i gruppi siano internamente il più omogenei possibili e il più possibile distanti fra loro, i.e. varianza tra le classi ampia e varianza entro piccola. In formule: \[\phi_*=argmax_a\frac{a^TBa}{a^tWa}\] dove \(B\) è la matrice di covarianza tra i gruppi e \(W\) quella entro.
Risolvendo questo problema di ottimizzazione si ricava: \[(W^{-1}B-\phi I_{G x G})a=0\]
Questo sistema di equazioni lineari omogenee ammette soluzioni non banali se e solo se \((W^{-1}B-\phi I)=0\), i.e. se \(\phi\) è un autovalore di \(W^{-1}B\) e \(a\) contiene il relativo autovettore. Algebricamente, questo implica che ci sono al massimo \(G-1\) soluzioni (rango di \(W^{-1}B\)) e, di conseguenza, \(G-1\) direzioni discriminanti. Ora, siccome \(\phi\) è proprio la quantità che vogliamo massimizzare, la funzione maggiormente discriminante sarà data dall’autovettore corrispondente all’autovalore maggiore di \(W^{-1}B\) e corrisponderà alla superficie che meglio separa un gruppo da tutti gli altri.
Nel caso di 2 classi avremo un solo autovalore non nullo e, quindi, una sola direzione discriminante. Per capire meglio il risultato di questo approccio possiamo infine rappresentare le osservazioni rispetto alla combinazione lineare \(a^Tx\), anche detta coordinata canonica, e sovraimporre il confine di decisione lineare.
Come si nota, le due classi non risultano molto separate. Probabilmente le discrete performance del modello sono dovute al fatto che le 2 categorie siano sbilanciate e che questa funzione discriminante permetta di classificare correttamente la maggior parte delle osservazioni di classe 1 in eccesso.verifica che le due formulazioni siano uguali
2.3. Quadratic Discriminant Analysis
Se facciamo cadere l’assunzione di omoschedasticità dei gruppi, l’approccio di analisi discriminante produce un funzione discriminante quadratica e, di conseguenza, anche la superficie di separazione non sarà più lineare.
risultati
2.4. Regressione Logistica
La regressione logistica è uno dei principali metodi utilizzati per problemi di classificazione che consente di effettuare la stima che l’evento interessato si realizzi (Y=1) per poi effettuare una classificazione delle osservazioni basata sul valore di probabilità stimata.
La stima della probabilità avviene sfruttando la funzione di ripartizione logistica tale che:
\[P(Y=1|X=x)=\frac{\exp{(\beta_0+\beta^Tx)}}{1+\exp{(\beta_0+\beta^Tx)}}\]
Applicando questo modello al nostro caso, i risultati ottenuti sono in linea con quelli osservati nei modelli precedenti anche se presentano delle performance leggermente inferiori avendo un’accuracy pari a 0.7388 e quindi un test error di 0.26117. Lo stesso discorso vale per i valori associati all’AUC che risulta leggermente inferiore agli altri con un valore di 0.8088. Osservando invece i valori di stima dei parametri, essi risultano tutti altamente significativi con unica eccezione per la variabile chlorides.
2.5. Regularised Methods
L’idea di questi metodi è quella di regolarizzare la stima dei parametri andando a “restringerli” verso lo zero. I metodi di regolarizzazione che abbiamo preso in considerazione per questo caso sono tre: ridge regression, logistic lasso ed elastic-net. Tutti e tre fanno uso del tuning parameter \(\lambda\) che influenza l’effetto della penalità, ma si differenziano per il diverso tipo di penalizzazione adottata. Nel seguito, l’addestramento del modello è condotto per cross-validation nel training set e i risultati di performance riportati sono calcolati poi sul test set.
La ridge regression utilizza una penalizzazione di tipo \(\ell_2\) per risolvere il problema di ottimizzazione dato da:
\[\min_{(\beta_0,\beta)\in\mathbb{R}^{p+1}}-{\biggl[ \frac{1}{n} \sum_{i=1}^{n}{y_i(\beta_0+x_i^T\beta)-\log{(1+\exp{(\beta_0+x_i^T\beta))}}} \biggr] + \lambda ||\beta||_2^2}\]
Questo tipo di formula permette di restringere i valori dei parametri verso zero senza mai annullarli del tutto. Al contrario il lasso, per valori sufficientemente alti di \(\lambda\)���, consente di rendere nulli i valori di alcuni (o tutti) i parametri svolgendo quindi anche funzione di selezione tramite l’utilizzo della penalizzazione di tipo \(\ell_1\):
\[\min_{(\beta_0,\beta)\in\mathbb{R}^{p+1}}-{\biggl[ \frac{1}{n} \sum_{i=1}^{n}{y_i(\beta_0+x_i^T\beta)-\log{(1+\exp{(\beta_0+x_i^T\beta))}}} \biggr] + \lambda ||\beta||_1}\]
Elastic-net, infine, costituisce un ibrido tra i primi due modelli che, introducendo una doppia penalità regolata dal parametro \(\alpha\), permette di combinare i vantaggi dei modelli precedenti:
\[\min_{(\beta_0,\beta)\in\mathbb{R}^{p+1}}-{\biggl[ \frac{1}{n} \sum_{i=1}^{n}{y_i(\beta_0+x_i^T\beta)-\log{(1+\exp{(\beta_0+x_i^T\beta))}}} \biggr] + \lambda[ \alpha ||\beta||_1 + (1-\alpha)||\beta||_2^2]}\]
Per tutti e tre i casi abbiamo utilizzato la 10-fold cross validation per individuare il valore ottimale di \(\lambda\) e quello di \(\alpha\) nel caso di elastic-net.
Analizzando i risultati ottenuti possiamo osservare come la ridge regression presenti il più alto valore di accuratezza e quindi inferiore test error rispetto agli altri due modelli e come tutti e tre abbiano performance inferiori rispetto alla regressione logistica. In tutti e tre i casi, infatti, i valori di \(\lambda\) ottenuti sono prossimi allo zero quasi a indicare che il miglior modello potrebbe essere quello non regolarizzato. Allo stesso modo la cross-validation ha individuato per elastic-net un valore di \(\alpha\) pari a 0.05 che indica che la principale penalizzazione utilizzata è quella associata alla ridge regression che in questo caso performa meglio del lasso. Per quanto riguarda l’accuratezza i risultati sono molto simili, con valori pari a, rispettivamente, 0.7305, 0.7345 e 0.7195 Questi risultati vengono confermati anche dalle ROC e soprattutto dalle AUC. I valori sono tutti molto vicini (0.80, 0.7995 e 0.7897) ma indicano comunque che la ridge regression ha un maggior potere predittivo.
Se osserviamo, però, i valori dei coefficienti, possiamo notare come la ridge regression per costruzione include tutti i parametri nell’analisi mentre il lasso permette di effettuare una selezione tra essi eliminando fixed.acidity, chlorides, density nonché type (vino bianco o rosso). Anche elastic-net, pur selezionato un valore di \(\alpha\) piuttosto piccolo, permette di effettuare una selezione dei parametri eliminando, però, solamente fixed.acidity e type. Nonostante le performance del lasso siano leggermente inferiori a quelle della ridge regression, potrebbe quindi valer la pena preferire il primo agli altri due in favore di un modello più semplice e di più facile interpretazione.
2.6. k Nearest Neighbor
Il K-Nearest Neighbour è un algoritmo di classificazione privo di parametri che si basa sul raggruppamento per classi in base al calcolo della distanza euclidea tra le osservazioni e sull’uso del “voto di maggioranza” per la classificazione. Nel caso preso in esame abbiamo fatto uso della 10-fold cross validation per la selezione del valore di K che definisce la dimensione del neighbourhood delle osservazioni. I valori sono stati tutti standardizzati in modo da evitare problemi dovuti a dimensioni di scala e la CV ci ha permesso di individuare un K ottimale pari a 23. I valori di accuracy sono in linea con gli altri finora ottenuti, tuttavia il KNN si posiziona al primo posto con un valore di 0.7572 e un’AUC di 0.8192.
2.7. Generalised Additive Models
I GAM sono modelli di tipo additivo che consentono di superare il vincolo di non linearità tramite l’introduzione di funzioni più flessibili nei regressori. Nel caso di classificazione binaria la forma assunta da un modello di tipo additivo è la seguente:
\[\log{ \biggl( \frac{P(Y=1|X)}{1-P(Y=1|X) } \biggr)}=\alpha + \sum_{k=1}^p{f_k(X)}\]
Nel caso preso in esame le funzioni utilizzate sono quelle di tipo smoothing spline e abbiamo fatto uso della 10-fold cross validation per selezionare il valore ottimale dei gradi di libertà utilizzando una ricerca a griglia con df da 1 a 4. Il valore selezionato è pari a 4. Possiamo osservare dai valori di accuratezza/test error che il modello così costruito presenta la miglior capacità predittiva rispetto agli altri finora analizzati con un’accuratezza dello 0.7658 confermata altresì dall’AUC di 0.8288.